Quality Control
jeff karnsomprot
2025-06-05
1 Objective
The goal of this session is to get a clean macula densa dataset to work with.
2 Load Packages
options(future.globals.maxSize = 74 * 1024^3) # 55 GB
getOption("future.globals.maxSize") #59055800320## [1] 794568949764 Subset out non MD cells.
SO2 <- subset(SO1, idents = c("4", "5", "6", "7"), invert = TRUE)
SO3 <- subset(SO2, features = grep("^mt-|^Gm|Rik|^Rp", rownames(SO2), invert = TRUE, value = TRUE))
DimPlot(SO2)
SO3 <- SCTransform(SO3) %>%
RunPCA() %>%
FindNeighbors(dims = 1:15) %>%
FindClusters(resolution = 2) %>%
RunUMAP(dims = 1:15)## Running SCTransform on assay: RNA## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 14658 by 11529## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells## There are 2 estimated thetas smaller than 1e-07 - will be set to 1e-07## Found 302 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 14658 genes## Computing corrected count matrix for 14658 genes## Calculating gene attributes## Wall clock passed: Time difference of 15.0815 secs## Determine variable features## Centering data matrix## Place corrected count matrix in counts slot## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Different cells and/or features from existing assay SCT## Set default assay to SCT## PC_ 1
## Positive: Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Prdx5, Foxq1, Klk1, Cox6c
## Wfdc15b, Krt7, Ly6a, Ckb, Slc25a5, Cox5b, Wfdc2, Atp5g1, Cldn19, Ggt1
## Ndufa4, Cox4i1, Atp1a1, Chchd10, Atp5b, Cox8a, Gadd45g, Cox6b1, Cox7a2, Uqcr11
## Negative: Pappa2, Zfand5, Aard, Robo2, Pde10a, Wwc2, Itga4, Nadk2, S100g, Nos1
## Ramp3, Neat1, Sgms2, Ptgs2, Col4a4, Irx1, Col4a3, Tmem158, Mir6236, Ranbp3l
## Itprid2, Bmp3, Cdkn1c, Camk2d, Dctd, Mcub, Sdc4, Etnk1, Peg3, Zbtb20
## PC_ 2
## Positive: Fos, Junb, Jun, Egr1, Btg2, Hspa1b, Hspa1a, Atf3, Zfp36, Ier2
## Fosb, Socs3, Klf6, Klf2, Jund, Dusp1, Dnajb1, Rhob, Gadd45g, Tsc22d1
## Gadd45b, Klf4, H3f3b, Ubc, Cebpd, Actb, H1f2, H2bc4, Ppp1r15a, Hspb1
## Negative: Pappa2, Egf, Slc12a1, CT010467.1, Sfrp1, Mir6236, Etnk1, Umod, Wnk1, Sec14l1
## Robo2, Hsp90b1, Pde10a, Car15, Atp1b1, Col4a4, Aard, Wwc2, Mal, Itprid2
## Mgst1, Aebp1, Col4a3, Bmp3, Trim2, Nos1, Sgms2, Zbtb20, Pth1r, Tfcp2l1
## PC_ 3
## Positive: Fth1, Ubb, Ftl1, Ldhb, Fxyd2, Car15, Prdx1, Mt1, Eif1, Cd63
## Cd9, Mpc2, Mgst1, Clu, Bsg, Tmem213, Aldoa, Mdh1, Ramp3, S100a1
## Selenow, Spp1, Tmem176b, Tspo, H3f3b, Wfdc2, Atpif1, Lgals3, Tmem59, Tmbim6
## Negative: Egf, Mir6236, CT010467.1, Umod, Tmem52b, Neat1, Nme7, Malat1, Slc12a1, Kcnq1ot1
## Lars2, Etnk1, Wnk1, Dst, Sult1d1, Foxq1, Slc5a3, Atrx, Atp1b1, Pnisr
## Fos, Syne2, Gls, Ivns1abp, Hnrnpa2b1, Srrm2, Chka, Col4a4, Paxbp1, Zbtb20
## PC_ 4
## Positive: Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Neat1, Ptger3, Ckb
## Fgf9, Egfl6, Ivns1abp, Defb1, Mfsd4a, Car4, Fkbp11, Chchd10, Plet1, Hspa1b
## Atp5md, Igfbp5, Cox6c, Tmem213, Atp5k, CT010467.1, Chka, Mgst3, Lpl, Atp1a1
## Negative: Pappa2, Aard, Tmem52b, Umod, Sult1d1, Tmem158, Egf, Foxq1, Mcub, Ptgs2
## Dctd, Wwc2, Iyd, Ramp3, S100g, Ptprz1, Cd9, Car15, Hsp90b1, Tmsb4x
## Defb42, Cdkn1c, Wnk1, Pth1r, Tagln2, Lcn2, Il1f6, Clu, Ctsc, Bmp2
## PC_ 5
## Positive: Hspa1b, Hspa1a, Umod, Egf, Jun, Cldn10, Pappa2, Car15, Klk1, Fos
## Sfrp1, Fth1, Ier3, Clcnkb, Atp1a1, Slc12a1, Cited2, Hspa8, Pik3r1, Cdkn1c
## Robo2, Kng2, Aard, Gsta4, Ptger3, Fau, Irx1, Ftl1, Uba52-ps, Tmem37
## Negative: S100g, Actb, Tmem52b, Igfbp7, Cryab, S100a6, Tacstd2, Tmsb10, Tmsb4x, Atf3
## Lcn2, Sult1d1, Foxq1, Crip1, Egr1, Ppia, Atp5md, Lgals1, Mir6236, Gadd45b
## Klf4, Ifitm3, Ccnd1, Ccn1, Hbegf, Serf2, Spp1, Nupr1, Cdkn1a, Krt7## Computing nearest neighbor graph## Computing SNN## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11529
## Number of edges: 359040
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7101
## Number of communities: 25
## Elapsed time: 1 seconds## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 13:27:36 UMAP embedding parameters a = 0.9922 b = 1.112## 13:27:36 Read 11529 rows and found 15 numeric columns## 13:27:36 Using Annoy for neighbor search, n_neighbors = 30## 13:27:36 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 13:27:37 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpIFjmQt/file54ff5c8bcfcc
## 13:27:37 Searching Annoy index using 1 thread, search_k = 3000
## 13:27:39 Annoy recall = 100%
## 13:27:39 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:27:40 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:27:40 Commencing optimization for 200 epochs, with 484894 positive edges
## 13:27:40 Using rng type: pcg
## 13:27:42 Optimization finished
markers.to.plot1 <- c("Lrp2", # PT
"Slc5a12", # PT-S1
"Slc13a3", # PT-S2
"Slc16a9", # PT-S3
"Havcr1", # Injured PT
"Epha7", # dTL
"Cryab", # dTL
"Cdh13", # dTL1
"Slc14a2", # dTL2
"Slc12a1", # TAL
"Umod", # TAL, DCT1
"Egf", # TAL, DCT1,
"Cldn10", # TAL
"Cldn16", # TAL
"Nos1", # MD
"Slc12a3", # DCT
"Pvalb", # DCT1
"Slc8a1", # DCT2, CNT
"Aqp2", # PC
"Slc4a1", # IC-A
"Slc26a4", # IC-B
"Nphs1", # Podo
"Ncam1", # PEC
"Flt1", # Endo
"Emcn", # Glom Endo
"Kdr", # Capillary Endo
"Pdgfrb", # Perivascular
"Pdgfra", # Fib
"Piezo2", # Mesangial
"Acta2", # Mural
"Ptprc", # Immune
"Cd74", # Macrophage
"Skap1", # B/T Cells
"Upk1b", # Uro
"Top2a" # Proliferation
)
DotPlot(SO3,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()## Warning: Could not find Slc5a12 in the default search locations, found in 'RNA'
## assay instead## Warning: Could not find Kdr in the default search locations, found in 'RNA'
## assay instead## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b
SO4 <- SCTransform(SO4) %>%
RunPCA() %>%
FindNeighbors(dims = 1:16) %>%
FindClusters(resolution = 2) %>%
RunUMAP(dims = 1:16)## Running SCTransform on assay: RNA## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 14577 by 11448## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells## There are 2 estimated thetas smaller than 1e-07 - will be set to 1e-07## Found 292 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 14577 genes## Computing corrected count matrix for 14577 genes## Calculating gene attributes## Wall clock passed: Time difference of 17.54207 secs## Determine variable features## Centering data matrix## Place corrected count matrix in counts slot## Warning: Different cells and/or features from existing assay SCT## Set default assay to SCT## PC_ 1
## Positive: Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Klk1, Prdx5, Cox6c, Foxq1
## Ly6a, Wfdc15b, Wfdc2, Slc25a5, Krt7, Cox5b, Ckb, Atp5g1, Cldn19, Cox4i1
## Ndufa4, Atp1a1, Ggt1, Cox8a, Atp5b, Chchd10, Cox6b1, Cox7a2, Uqcrq, Uqcr11
## Negative: Pappa2, Zfand5, Aard, Robo2, Wwc2, Pde10a, Mir6236, Itga4, Nadk2, Neat1
## Nos1, S100g, Col4a4, Ramp3, Sgms2, Col4a3, Ptgs2, Irx1, Tmem158, Itprid2
## Ranbp3l, Bmp3, Camk2d, Sdc4, Etnk1, Cdkn1c, Zbtb20, Srrm2, Hnrnpa2b1, Peg3
## PC_ 2
## Positive: Fos, Junb, Jun, Hspa1b, Hspa1a, Btg2, Egr1, Zfp36, Atf3, Ier2
## Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Tsc22d1, Rhob, Gadd45g
## H3f3b, Gadd45b, Ubc, Cebpd, Mt1, Mt2, Actb, H2bc4, H1f2, Klf4
## Negative: Mir6236, CT010467.1, Egf, Pappa2, Slc12a1, Umod, Etnk1, Wnk1, Nme7, Atp1b1
## Sfrp1, Robo2, Sptbn1, Col4a4, Sec14l1, Hsp90b1, Zbtb20, Dst, Itprid2, Tmem52b
## Lars2, Pde10a, Kcnq1ot1, App, Mal, Atrx, Col4a3, Tfcp2l1, Syne2, Nfe2l1
## PC_ 3
## Positive: Egf, Mir6236, CT010467.1, Umod, Tmem52b, Neat1, Nme7, Fos, Jun, Malat1
## Slc12a1, Junb, Kcnq1ot1, Lars2, Egr1, Sult1d1, Foxq1, Etnk1, Dst, Wnk1
## Slc5a3, Atp1b1, Atf3, Fosb, Zfp36, Btg2, Atrx, Klf6, Pnisr, Ivns1abp
## Negative: Fth1, Ubb, Ftl1, Fxyd2, Ldhb, Car15, Prdx1, Cd63, Cd9, Eif1
## Mpc2, Mgst1, Mt1, Tmem213, Bsg, Clu, Mdh1, Aldoa, Itm2b, Ramp3
## Selenow, Tmem59, S100a1, Spp1, Tmem176b, Tspo, Tmbim6, Tle5, Atpif1, Wfdc2
## PC_ 4
## Positive: Pappa2, Aard, Tmem52b, Umod, Egf, Sult1d1, Tmem158, Foxq1, Mcub, Ptgs2
## Dctd, Wwc2, Iyd, Car15, Ramp3, Ptprz1, Cd9, Hsp90b1, S100g, Wnk1
## Pth1r, Cdkn1c, Defb42, Clu, Ctsc, Tmsb4x, Tagln2, Bmp2, Col4a3, Tsc22d1
## Negative: Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Ptger3, Neat1, Ckb
## Fgf9, Egfl6, Ivns1abp, Mfsd4a, Defb1, Car4, CT010467.1, Plet1, Fkbp11, Atp5md
## Chchd10, Cox6c, Atp5k, Igfbp5, Tmem213, Mgst3, Chka, Atp5mpl, Avpr1a, Lpl
## PC_ 5
## Positive: CT010467.1, Hspa1b, Hspa1a, Klk1, Pappa2, Car15, Cldn10, Mir6236, Fth1, Lars2
## Jun, Hspa8, Wfdc2, Itm2b, Ftl1, Fau, Uba52-ps, Ly6a, Sfrp1, Aard
## Pik3r1, Mal, Eef1a1, Atp1a1, Ptger3, Atp1b1, Tspo, Hspb1, Tmem176a, Gpx4
## Negative: S100g, Actb, Tmem52b, Sdc4, Abhd2, Egr1, Uroc1, Tmsb10, Ppia, Atf3
## Serf2, Slc39a1, Alkbh5, Ndufb1-ps, Cebpd, Igfbp7, Cox6c, Atp5md, Ramp3, Ndufa3
## Wfdc15b, Kdm6b, Atp5k, Atp5e, Rbm47, Ldhb-ps, Ranbp3l, Ppm1h, Atp5mpl, Rcan1## Computing nearest neighbor graph## Computing SNN## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11448
## Number of edges: 356414
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7087
## Number of communities: 26
## Elapsed time: 0 seconds## 13:28:13 UMAP embedding parameters a = 0.9922 b = 1.112## 13:28:13 Read 11448 rows and found 16 numeric columns## 13:28:13 Using Annoy for neighbor search, n_neighbors = 30## 13:28:13 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 13:28:13 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpIFjmQt/file54ff32e1dca1
## 13:28:13 Searching Annoy index using 1 thread, search_k = 3000
## 13:28:15 Annoy recall = 100%
## 13:28:15 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:28:16 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:28:16 Commencing optimization for 200 epochs, with 483954 positive edges
## 13:28:16 Using rng type: pcg
## 13:28:18 Optimization finished
SO4@meta.data$seurat_clusters <- factor(
SO4@meta.data$seurat_clusters,
levels = c("0", "1", "2", "3", "4", "6", "8", "9", "10", "12", "15", "22", "18",
"19", "21", "16", "13", "20", "23", "17", "14", "5", "7", "11", "24", "25")
)
Idents(SO4) <- SO4$seurat_clusters
## Simplifying into a few types
markers.to.plot2 <- c("Pappa2",
"Nos1",
"Egf",
"Umod",
"Fos",
"Junb",
"Cxcl10",
"Isg15",
"Sult1d1",
"Cldn19",
"Foxq1"
)
DimPlot(SO4)
DotPlot(SO4,
features = markers.to.plot2,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "seurat_clusters",
col.max = 2.5)+
coord_flip()
Idents(SO4) <- "seurat_clusters"
all_markers <- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 0.5)## Calculating cluster 0## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session## Calculating cluster 1## Warning in FindMarkers.default(object = data.use, cells.1 = cells.1, cells.2 =
## cells.2, : No features pass logfc.threshold threshold; returning empty
## data.frame## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 6## Calculating cluster 8## Calculating cluster 9## Calculating cluster 10## Calculating cluster 12## Calculating cluster 15## Calculating cluster 22## Calculating cluster 18## Calculating cluster 19## Calculating cluster 21## Calculating cluster 16## Calculating cluster 13## Calculating cluster 20## Calculating cluster 23## Calculating cluster 17## Calculating cluster 14## Calculating cluster 5## Calculating cluster 7## Calculating cluster 11## Calculating cluster 24## Calculating cluster 25all_markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
slice_head(n = 3) %>%
ungroup() -> top3
DoHeatmap(SO4, features = top3$gene) + NoLegend()
5 Assigning Subclasses for clusters
SO4@meta.data <- SO4@meta.data %>%
mutate(subclass_MD = dplyr::case_when(
seurat_clusters %in% c(0,1,2,4,6,10,21,18,19,12,22) ~ "type_1a",
seurat_clusters %in% c(3,9,15,21,8) ~ "type_1b",
seurat_clusters %in% c(7,24,25,5,11) ~ "type_2a",
seurat_clusters %in% c(14,17) ~ "type_2b",
seurat_clusters %in% c(16,13) ~ "type_3a",
seurat_clusters %in% c(20) ~ "type_3b",
seurat_clusters == 23 ~ "type_4",
))
SO4@meta.data$subclass_MD <- factor(
SO4@meta.data$subclass_MD,
levels = c("type_1a","type_1b", "type_2a", "type_2b","type_3a","type_3b","type_4")
)
SO4@meta.data <- SO4@meta.data %>%
mutate(subclass2_MD = dplyr::case_when(
seurat_clusters %in% c(0,1,2,3,4,6,8,9,10,21,15,22,18,19,21,12) ~ "type_1",
seurat_clusters %in% c(7,24,25,5,11,14,17) ~ "type_2",
seurat_clusters %in% c(16,13,20) ~ "type_3",
seurat_clusters %in% c(23) ~ "type_4"
))
SO4@meta.data$subclass2_MD <- factor(
SO4@meta.data$subclass2_MD,
levels = c("type_1", "type_2","type_3","type_4")
)
Idents(SO4) <- SO4@meta.data$subclass_MD
DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)
library(plotly)
dim_plot <- DimPlot(SO4, reduction = "umap",group.by = "seurat_clusters") +
ggtitle("UMAP Dimensional Reduction")
interactive_plot <- ggplotly(dim_plot)
# Display the interactive plot
interactive_plot# Assign your markers as a vector
genes_to_plot <- c(
"Fos", "Socs3", "Gadd45b", "Wee1", "Hspa1a", "Sat1", # TYPE 3
"Cxcl12", "Itpr2", "Bmp4", "Casr", "Grin2c", "Irx3", "Rap1gap", "App", "Wwc1", # TYPE 2
"Syt5", "Syn3", "Cacna1d", "Slc6a7", "Robo2", "Begain", # TYPE 1
"Pappa1", "Nos1", "Ptgs2", "Bmp2", "Atp2a3" # SUBCLUSTERS within type 1
)
DoHeatmap(SO4, features = genes_to_plot, group.by = "seurat_clusters")## Warning in DoHeatmap(SO4, features = genes_to_plot, group.by =
## "seurat_clusters"): The following features were omitted as they were not found
## in the scale.data slot for the SCT assay: Pappa1, Wwc1, Irx3, Grin2c
6 Viewing genes that differentiate each Type
markers.to.plot3 <- c("Pappa2",
"Egf",
"Umod",
"Fos",
"Junb",
"Cxcl10",
"Isg15"
)
DotPlot(SO4,
features = markers.to.plot3,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "subclass2_MD",
col.max = 2.5)+
coord_flip()
DotPlot(SO4,
features = markers.to.plot2,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "subclass_MD",
col.max = 2.5)+
coord_flip()
Idents(SO4) <- "seurat_clusters"
all_markers <- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 0.5)## Calculating cluster 0## Calculating cluster 1## Warning in FindMarkers.default(object = data.use, cells.1 = cells.1, cells.2 =
## cells.2, : No features pass logfc.threshold threshold; returning empty
## data.frame## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 6## Calculating cluster 8## Calculating cluster 9## Calculating cluster 10## Calculating cluster 12## Calculating cluster 15## Calculating cluster 22## Calculating cluster 18## Calculating cluster 19## Calculating cluster 21## Calculating cluster 16## Calculating cluster 13## Calculating cluster 20## Calculating cluster 23## Calculating cluster 17## Calculating cluster 14## Calculating cluster 5## Calculating cluster 7## Calculating cluster 11## Calculating cluster 24## Calculating cluster 25# Order by avg_log2FC
all_markers <- all_markers %>%
arrange(desc(avg_log2FC)) %>%
select(gene, everything())
# Split by cluster (ident column)
marker_list <- split(all_markers, all_markers$cluster)
# Sort names alphanumerically
sorted_cluster_names <- names(marker_list)[order(as.numeric(names(marker_list)))]
# Create a workbook
wb <- createWorkbook()
# Add each cluster as a new worksheet
for (cluster_name in sorted_cluster_names) {
addWorksheet(wb, sheetName = cluster_name)
writeData(wb, sheet = cluster_name, marker_list[[cluster_name]])
}
date <- format(Sys.Date(), "%Y%m%d")
# Save workbook
saveWorkbook(wb, here("jk_code", paste0(date, "_", "_SO4_FindAllMarkers_By_Cluster.xlsx")), overwrite = TRUE)
# Saving Seurat Object
save(SO4, file = here("jk_code", "JK_clean_MD.rds"))all_markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
slice_head(n = 3) %>%
ungroup() -> top3
DoHeatmap(SO4, features = top3$gene) + NoLegend()## Warning in DoHeatmap(SO4, features = top3$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Nubp2